{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Forecasting with sktime - appendix: forecasting, supervised regression, and pitfalls in confusing the two\n",
    "\n",
    "This notebook provides some supplementary explanation about the relation between forecasting as implemented in `sktime`, and the very common supervised prediction tasks as supported by `scikit-learn` and similar toolboxes.\n",
    "\n",
    "Key points discussed in this notebook:\n",
    "\n",
    "* forecasting is not the same as supervised prediction\n",
    "* even though forecasting can be \"solved\" by algorithms for supervised prediction, this is indirect and requires careful composition\n",
    "* from an interface perspective, this is correctly formulated as \"reduction\", i.e., use of a supervised predictor as a component within a forecaster\n",
    "* there are a number of pitfalls if this is manually done - such as, over-optimistic performance evaluation, information leakage, or \"predicting the past\" type errors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# general imports\n",
    "import numpy as np\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The pitfalls of mis-diagnosing forecasting as supervised regression\n",
    "\n",
    "A common mistake is to mis-identify a forecasting problem as supervised regression - after all, in both we predict numbers, so surely this must be the same thing?\n",
    "\n",
    "Indeed we predict numbers in both, but the set-up is different:\n",
    "\n",
    "* in supervised regression, we predict *label/target variables* from *feature variables*, in a cross-sectional set-up. This is after training on label/feature examples.\n",
    "* in forecasting, we predict *future values* from *past values*, of *the same variable*, in a temporal/sequential set-up. This is after training on the past.\n",
    "\n",
    "In the common data frame representation:\n",
    "\n",
    "* in supervised regression, we predict entries in a column from other columns. For this, we mainly make use of the statistical relation between those columns, leart from examples of complete rows. The rows are all assumed exchangeable.\n",
    "* in forecasting, we predict new rows, assuming temporal ordering in the rows. For this, we mainly make use of the statistical relation between previous and subsequent rows, learnt from the example of the observed sequence of rows. The rows are not exchangeable, but in temporal sequence.\n",
    "\n",
    "TODO: add a nice picture on what is predicted from what, arrows and all. Contributions are welcome."
   ]
  },
  {
   "source": [
    "### Pitfall 1: over-optimism in performance evaluation, false confidence in \"broken\" forecasters\n",
    "\n",
    "Confusing the two tasks may lead to information leakage, and over-optimistic performance evaluation. This is because in supervised regression the ordering of rows does not matter, and train/test split is usually performed uniformly. In forecasting, the ordering does matter, both in training and in evaluation.\n",
    "\n",
    "As subtle as it seems, this may have major practical consequences - since it can lead to the mistaken belief that a \"broken\" method is performant, which can cause damage to health, property, and other assets in real-life deployment.\n",
    "\n",
    "The example below shows \"problematic\" performance estimation, when mistakenly using the regression evaluation workflow for forecasting."
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-10T16:07:06.641893Z",
     "iopub.status.busy": "2021-04-10T16:07:06.618023Z",
     "iopub.status.idle": "2021-04-10T16:07:06.787235Z",
     "shell.execute_reply": "2021-04-10T16:07:06.787715Z"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "from sktime.datasets import load_airline\n",
    "from sktime.forecasting.model_selection import temporal_train_test_split\n",
    "from sktime.utils.plotting import plot_series"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = load_airline()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_train, y_test = train_test_split(y)\n",
    "plot_series(y_train.sort_index(), y_test.sort_index(), labels=[\"y_train\", \"y_test\"]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This leads to leakage:\n",
    "\n",
    "> The data you are using to train a machine learning algorithm happens to have the information you are trying to predict.\n",
    "\n",
    "But `train_test_split(y, shuffle=False)` works, which is what `temporal_train_test_split(y)` does in `sktime`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-10T16:07:06.861616Z",
     "iopub.status.busy": "2021-04-10T16:07:06.861025Z",
     "iopub.status.idle": "2021-04-10T16:07:07.016695Z",
     "shell.execute_reply": "2021-04-10T16:07:07.017191Z"
    }
   },
   "outputs": [],
   "source": [
    "y_train, y_test = temporal_train_test_split(y)\n",
    "plot_series(y_train, y_test, labels=[\"y_train\", \"y_test\"]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pitfall 2: obscure data manipulations, brittle boilerplate code to apply regressors\n",
    "\n",
    "It is common practice to apply supervised regressors after transforming the data for forecasting, through lagging - for example, in auto-regressive reduction strategies.\n",
    "\n",
    "Two important pitfalls appear right at the start:\n",
    "\n",
    "* a lot of boilerplate code has to be written to transform the data to make it ready for fitting - this is highly error prone\n",
    "* there are a number of implicit hyper-parameters here, such as window and lag size. If done without caution, these are not explicit or tracked in the experiment, which can lead to \"p-value hacking\".\n",
    "\n",
    "Below is an example of such boilerplate code to demonstrate this. The code is closely modelled on the R code used in the [M4 competition](https://github.com/Mcompetitions/M4-methods):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# suppose we want to predict 3 years ahead\n",
    "fh = np.arange(1, 37)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-10T16:07:07.023411Z",
     "iopub.status.busy": "2021-04-10T16:07:07.022907Z",
     "iopub.status.idle": "2021-04-10T16:07:07.024499Z",
     "shell.execute_reply": "2021-04-10T16:07:07.024998Z"
    }
   },
   "outputs": [],
   "source": [
    "# slightly modified code from the M4 competition\n",
    "def split_into_train_test(data, in_num, fh):\n",
    "    \"\"\"\n",
    "    Splits the series into train and test sets.\n",
    "\n",
    "    Each step takes multiple points as inputs\n",
    "    :param data: an individual TS\n",
    "    :param fh: number of out of sample points\n",
    "    :param in_num: number of input points for the forecast\n",
    "    :return:\n",
    "    \"\"\"\n",
    "    train, test = data[:-fh], data[-(fh + in_num) :]\n",
    "    x_train, y_train = train[:-1], np.roll(train, -in_num)[:-in_num]\n",
    "    x_test, y_test = test[:-1], np.roll(test, -in_num)[:-in_num]\n",
    "    #     x_test, y_test = train[-in_num:], np.roll(test, -in_num)[:-in_num]\n",
    "\n",
    "    # reshape input to be [samples, time steps, features]\n",
    "    # (N-NF samples, 1 time step, 1 feature)\n",
    "    x_train = np.reshape(x_train, (-1, 1))\n",
    "    x_test = np.reshape(x_test, (-1, 1))\n",
    "    temp_test = np.roll(x_test, -1)\n",
    "    temp_train = np.roll(x_train, -1)\n",
    "    for x in range(1, in_num):\n",
    "        x_train = np.concatenate((x_train[:-1], temp_train[:-1]), 1)\n",
    "        x_test = np.concatenate((x_test[:-1], temp_test[:-1]), 1)\n",
    "        temp_test = np.roll(temp_test, -1)[:-1]\n",
    "        temp_train = np.roll(temp_train, -1)[:-1]\n",
    "\n",
    "    return x_train, y_train, x_test, y_test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-10T16:07:07.028644Z",
     "iopub.status.busy": "2021-04-10T16:07:07.028108Z",
     "iopub.status.idle": "2021-04-10T16:07:07.029820Z",
     "shell.execute_reply": "2021-04-10T16:07:07.030335Z"
    }
   },
   "outputs": [],
   "source": [
    "# here we split the time index, rather than the actual values,\n",
    "# to show how we split the windows\n",
    "feature_window, target_window, _, _ = split_into_train_test(\n",
    "    np.arange(len(y)), 10, len(fh)\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To better understand the prior data transformation, we can look at how we can split the training series into windows. Here we show the generated windows expressed as integer indices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-10T16:07:07.033532Z",
     "iopub.status.busy": "2021-04-10T16:07:07.033064Z",
     "iopub.status.idle": "2021-04-10T16:07:07.035278Z",
     "shell.execute_reply": "2021-04-10T16:07:07.035836Z"
    }
   },
   "outputs": [],
   "source": [
    "feature_window[:5, :]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-10T16:07:07.038874Z",
     "iopub.status.busy": "2021-04-10T16:07:07.038346Z",
     "iopub.status.idle": "2021-04-10T16:07:07.040421Z",
     "shell.execute_reply": "2021-04-10T16:07:07.040933Z"
    }
   },
   "outputs": [],
   "source": [
    "target_window[:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-10T16:07:07.044178Z",
     "iopub.status.busy": "2021-04-10T16:07:07.043592Z",
     "iopub.status.idle": "2021-04-10T16:07:07.046140Z",
     "shell.execute_reply": "2021-04-10T16:07:07.046641Z"
    }
   },
   "outputs": [],
   "source": [
    "# now we can split the actual values of the time series\n",
    "x_train, y_train, x_test, y_test = split_into_train_test(y.values, 10, len(fh))\n",
    "print(x_train.shape, y_train.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-10T16:07:07.083439Z",
     "iopub.status.busy": "2021-04-10T16:07:07.082648Z",
     "iopub.status.idle": "2021-04-10T16:07:07.167067Z",
     "shell.execute_reply": "2021-04-10T16:07:07.167558Z"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.ensemble import RandomForestRegressor\n",
    "\n",
    "model = RandomForestRegressor()\n",
    "model.fit(x_train, y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To reiterate the potential pitfalls here:\n",
    "\n",
    "> The manual requires a lot of hand-written code which is often error-prone, not modular and not tuneable.\n",
    "\n",
    "> These steps involve a number of implicit hyper-parameters:\n",
    "> * the way you slice the time series into windows (e.g. the window length)\n",
    "> * the way you generate forecasts (recursive strategy, direct strategy, other hybrid strategies)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pitfall 3: Given a fitted regression algorithm, how can we generate forecasts?\n",
    "\n",
    "The next important pitfall comes at the end:\n",
    "\n",
    "if making predictions along the \"manual route\" for supervised regressors, the supervised regressor's outputs have to be transformed back into forecasts. This is easily forgotten, and invites errors in forecasts and evaluation (see pitfall no.1) - especially, if one does not cleanly keep track of which data is known at what time, or how to invert the transformation made in fitting.\n",
    "\n",
    "A naive user might now proceed like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-10T16:07:07.171836Z",
     "iopub.status.busy": "2021-04-10T16:07:07.171233Z",
     "iopub.status.idle": "2021-04-10T16:07:07.173018Z",
     "shell.execute_reply": "2021-04-10T16:07:07.173506Z"
    }
   },
   "outputs": [],
   "source": [
    "print(x_test.shape, y_test.shape)\n",
    "\n",
    "# add back time index to y_test\n",
    "y_test = pd.Series(y_test, index=y.index[-len(fh) :])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-10T16:07:07.180067Z",
     "iopub.status.busy": "2021-04-10T16:07:07.177212Z",
     "iopub.status.idle": "2021-04-10T16:07:07.186432Z",
     "shell.execute_reply": "2021-04-10T16:07:07.186959Z"
    }
   },
   "outputs": [],
   "source": [
    "y_pred = model.predict(x_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sktime.performance_metrics.forecasting import mean_absolute_percentage_error\n",
    "\n",
    "mean_absolute_percentage_error(pd.Series(y_pred, index=y_test.index), y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So easy, so wrong ... but what's the problem here? It's a bit subtle and not easy to spot:\n",
    "\n",
    "> We actually don't make a multi-step-ahead forecast up to the 36th step ahead. Instead, we make 36 single-step-ahead forecasts always using the most recent data. But that's a solution to a different learning task!\n",
    "\n",
    "To fix this problem, we could write some code to do this recursively as in the M4 competition:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-10T16:07:07.228733Z",
     "iopub.status.busy": "2021-04-10T16:07:07.219457Z",
     "iopub.status.idle": "2021-04-10T16:07:07.421143Z",
     "shell.execute_reply": "2021-04-10T16:07:07.421651Z"
    }
   },
   "outputs": [],
   "source": [
    "# slightly modified code from the M4 study\n",
    "predictions = []\n",
    "last_window = x_train[-1, :].reshape(1, -1)  # make it into 2d array\n",
    "\n",
    "last_prediction = model.predict(last_window)[0]  # take value from array\n",
    "\n",
    "for i in range(len(fh)):\n",
    "    # append prediction\n",
    "    predictions.append(last_prediction)\n",
    "\n",
    "    # update last window using previously predicted value\n",
    "    last_window[0] = np.roll(last_window[0], -1)\n",
    "    last_window[0, (len(last_window[0]) - 1)] = last_prediction\n",
    "\n",
    "    # predict next step ahead\n",
    "    last_prediction = model.predict(last_window)[0]\n",
    "\n",
    "y_pred_rec = pd.Series(predictions, index=y_test.index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sktime.performance_metrics.forecasting import mean_absolute_percentage_error\n",
    "\n",
    "mean_absolute_percentage_error(pd.Series(y_pred, index=y_test.index), y_test)"
   ]
  },
  {
   "source": [
    "To summarize the potential pitfalls here:\n",
    "\n",
    "> Obtaining regressor predictions and converting them back into forecasts is non-trivial and error prone:\n",
    "> * some boilerplate code needs to be written, which just as in pitfall no.2 introduces potential for problems\n",
    "> * it isn't exactly obvious that this boilerplate code had to be written in the first place, creating a subtle failure point"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "source": [
    "### How does `sktime` help avoid the above pitfalls?\n",
    "\n",
    "`sktime` mitigates the above pitfalls by:\n",
    "\n",
    "* the unified interface for forecasters - any strategy to produce forecasts is a forecaster. Through the unified interface, forecasters are directly compatible with deployment and evaluation workflows appropriate for forecasters.\n",
    "* its declarative specification interface that minimizes boilerplate code - it's minimized to the bare necessities to tell `sktime` which forecaster you want to build\n",
    "\n",
    "Nevertheless, `sktime` aims to be flexible, and tries to avoid to railroad the user into specific methodological choices."
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.neighbors import KNeighborsRegressor\n",
    "\n",
    "from sktime.forecasting.compose import make_reduction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# declarative forecaster specification - just two lines!\n",
    "regressor = KNeighborsRegressor(n_neighbors=1)\n",
    "forecaster = make_reduction(regressor, window_length=15, strategy=\"recursive\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "forecaster.fit(y)\n",
    "y_pred = forecaster.predict(fh)"
   ]
  },
  {
   "source": [
    "... and that's it!\n",
    "\n",
    "Note that there is no `x_train` or other boilerplate artefacts, since construction of the lagged features and other boilerplate code are taken care of by the forecaster internally.\n",
    "\n",
    "For more details on the `sktime` composition interface, refer to Section 3 of the main forecasting tutorial."
   ],
   "cell_type": "markdown",
   "metadata": {}
  }
 ],
 "metadata": {
  "hide_input": false,
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.10"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
